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1  INTRODUCTION 

Mechanical  engineers  are  faced  with  the  task  to  design  tools  which  will  operate  safely 
under  certain  mechanical  conditions,  in  certain  environments  and  for  certain  periods 
of  time.  The  cost  of  fabrication  and  design  has  to  be  low  and  the  construction  has 
to  meet  specific  demands  on  weight  and  space  requirements  etc. 

Before  manufacturing  the  tool,  the  designer  must  be  able  to  predict  its  behaviour. 
This  piediction  is  based  on  a  formulation  of  a  mathematical  model,  its  computational 
analysis,  experiments,  and  experience  with  existing  constructions  and  their  failures. 
Because  of  various  uncertainties  which  necessarily  occur,  the  goals  of  the  advanced 
design  analysis  (in  nuclear  industry,  aircraft  industry  etc.)  are  often  stipulated  in 
the  design  codes  (which  are  changing  over  time)  and  are,  at  least  by  parts,  company 
oriented.  The  question  of  the  principles  of  the  safety  is  directly  related  to  these 
codes. 

For  example,  in  the  design  code  [1]  used  in  military  aircraft  design,  it  is  required 
that  components  based  on  the  principles  of  ’'non-inspectable  slow  crack  growth" 
must  be  designed  under  the  assumptions  that. 


a)  the  as-fabricated  structure  contains  flaws  of  a  size  just  smaller  than  the  non¬ 
destructive  maximum  undetectable  flaw  size 

b)  the  flaws  are  assumed  to  exist  in  form  of  crack-like  defects  with  most  un¬ 
favourable  location  and  orientation 


The  design  code  requires  that  the  mathematical  formulation  and  its  computa¬ 
tional  analysis  must  reliably  and  conservatively  predict  both  the  sizes  of  the  growing 
cracks  and  the  residual  strength  of  the  component. 

Of  course  an  essential  part  of  these  principles  is  nonmathematical  as  for  example 
the  mechanical  principles  of  fatigue  crack  growth  (today  mostly  based  on  the  phe- 
nomencdogical  Paris  law  or  similar  laws).  For  various  aspects  we  refer  the  reader  to 
[2], 

There  obviously  must  be  a  balance  between  what  can  be  accurately  analyzed 
and  what  the  design  code  stipulates.  This  is  influenced  by  the  steady  progress 
of  computer  hardware  as  well  as  progress  in  the  numerical  methods.  Today  and 
still  in  the  future  more  and  more  complex  mathematical  problems  can  and  will  be 
understood  and  reliably  solved.  The  role  of  the  effectiveness  of  the  computations, 
although  very  important,  has  to  be  assessed  from  this  viewpoint,  especially  when 
the  effectiveness  could  drastically  depend  on  the  computer  architecture. 

The  role  of  numerical  analysis  has  to  be  seen  also  in  the  light  of  testing  the 
mathematical  models  (the  formulation)  by  comparing  the  computed  results  with 


6 


F  FA  TNI  992- 1 7 


experiments.  Here  it  is  essential  that  the  numerical  results  are  so  close  to  the 
exact  solution  of  the  mathematical  problem  that  any  observed  discrepancy  with 
experimental  results  is  only  due  to  the  mathematical  formulation  and  not  the  errors 
of  the  numerical  analysis.  This  means  that  a  practical  accurate  and  reliable  a- 
posteriori  error  estimate  for  any  data  of  interest  (i.e.  not  only  the  usual  energy 
norm  estimate)  is  imperative. 

The  FENcode  has  to  be  adaptive  in  the  sense  that  given  a  tolerance  for  the  error  in 
stress,  stress  intensity  factor  etc.  the  FENcode  automatically  increases  the  accuracy 
until  the  error  becomes  lower  than  requested  by  the  user.  The  increase  in  accuracy 
of  a  FE-solution  may  be  obtained  in  basically  three  different  ways. 


a)  The  /i-version.  Here  the  degrees  of  the  elements  are  fixed  and  the  program 
adaptively  refines  the  mesh  (uniformly  or  non-uniformly)  until  the  desired 
accuracy  is  achieved.  There  are  several  h-adaptive  codes.  We  refer  to  [3], [4] 
as  examples 

b)  The  p-version.  Here  the  mesh  is  fixed  and  the  element  degrees  p  are  adaptively 
increased  (possibly  non-uniformly)  until  the  desired  accuracy  is  achieved.  The 
programs  MSC/PROBE  (McN’eal  Schwendler  Corp.).  STRIPE  (The  Aeronau¬ 
tical  Research  Institute  of  Sweden),  Applied  Structure  (Rasna  Corp.)  belong 
to  this  category.  Although  in  principle,  the  method  gives  the  result  needed  for 
an\'  mesh,  practically  the  mesh  is  of  utmost  importance. 

c)  The  h  —  p  version.  Here  the  code  adaptively  refines  the  mesh  and  selects 
the  degrees  p.  The  code  PHLEX  (Computational  Mechanics  Corp..  .Austin) 
utilizes  the  h  —  p  version  of  FEM 


If  the  p-version  is  used  on  a  properl\  designed  mesh  it  essentially  produces  the 
accuracy  of  the  h  —  p  version.  By  following  simple  rules  for  mesh  design  one  almost 
achieves  the  exponential  convergence  rate  for  accuracies  of  interest.  In  what  will 
follow,  we  concentrate  on  the  p-version  of  FE.M  which  combined  with  properly  graded 
meshes  gives  a  convergence  rate  similar  to  that  of  the  h  —  p  version. 

One  of  the  major  advantages  of  the  p-version  is  that  in  sequential  computations 
for  increasing  p  a  sequence  of  solutions  for  the  parameter  of  interest  is  obtained 
from  which  one  can  make  very  reliable  assessments  of  the  accuracy  of  the  solution. 
In  principle  this  can  be  done  for  the  h-vevsion  too,  although  for  practical  three- 
dimensional  problems  the  number  of  degrees  of  freedom  quickly  become  prohibitively 
high. 

In  this  and  forthcoming  papers  we  will  especially  concentrate  on  problems  which 
are  relevant  to  the  criteria  of  fatigue  crack  initiation  and  stable  crack  growth  of  flaws 
in  metallic  structures. 
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The  mathematical  model  considered  is  the  Navier-Lame  equations  of  three- 
dimensional  linear  elasticity.  We  will  address  the  mathematical  models  and  the 
numerical  procedures  and  give  practical  examples  of  the  achievements.  W'e  will  con¬ 
centrate  here  on  the  cases  where  the  mathematical  model  is  given  in  a  deterministic 
way.  To  make  our  ideas  as  clear  as  possible  we  will  analyze  an  academic  example 
and  a  complex  fuselage  part. 

In  section  2  we  will  discuss  a  three-dimensional  elasticity  model  problem.  We 
concentrate  on  its  solution  in  a  neighbourhood  of  the  edges  and  vertices.  The  singu¬ 
lar  behaviour  of  the  solution  in  these  neighbourhoods  is  characterized  by  edge  stress 
intensity  functions  and  vertex  stress  intensity  factors.  We  show  that  the  p-version 
on  a  proper  mesh  together  with  use  of  extraction  procedures  leads  to  accurate  com¬ 
putations  of  these  functions  and  factors  as  well  as  stresses  with  virtually  guaranteed 
a-posteriori  error  estimation.  We  also  show  that  the  stress  intensity  functions  and 
factors  converge  with  the  same  rate  as  the  energy  error  of  the  finite  element  solution 
(i.e.  as  the  square  of  the  energy  norm)  and  that  the  convergence  is  exponential  in 
the  accuracy  range  of  practical  interest. 

In  the  section  3  we  adress  the  problem  of  reliable  stress  calculation  in  complex 
three-dimensional  components.  A  complex  fuselage  frame  is  analyzed  with  the  fa¬ 
tigue  design  problem  in  mind.  An  effective  method  for  the  design  of  the  radices  at 
fillets  is  described.  We  concentrate  on  the  reliability  of  the  computed  data  and  the 
a-posteriort  error  estimation  in  the  context  of  the  analysis  of  the  fuselage  frame  prob¬ 
lem.  Section  4  discusses  the  damage  tolerance  analysis  o{  complex  three-dimensional 
components.  The  fuselage  component  frame  with  cracks  of  different  sizes  is  analyzed. 
Section  5  adresses  the  computational  zispects  of  the  h  —  p  version.  It  describes  the 
essential  parts  of  the  computations  and  presents  basic  information  about  Cpu-times 
and  memory  needed.  The  sizes  of  the  problems  are  60000  to  1260000  degrees  of 
freedom. 

Section  6  summarices  the  basic  properties  of  the  h  —  p  version,  its  relibably  error 
estimation  and  practical  experience.  All  data  and  experience  is  based  on  the  use  of 
the  program  STRIPE  developed  at  the  Aeronautical  Research  Institute  of  Sweden. 

The  computational  methods  presented  here  are  very  well  covered  by  a  math¬ 
ematical  theory.  In  the  following  papers  [-o],  [6],[7],[8].[9].[10]  we  will  elaborate  in 
detail  on  these  methods  and  their  computational  aspects. 
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2  A  MODEL  PROBLEM  AND  ITS  ANALYSIS 


In  order  to  exemplify  the  character  of  the  solutions  to  problems  in  three-dimensional 
elastomechanics,  we  consider  the  folded  plate  £l  shown  in  figure  1.  The  plate  is  loaded 
with  parabolic  sL .  ,r  in  the  symmetry  section. 


Figure  1:  Folded  plate  loaded  with  parabolic  shear  in  symmetry  section.  Plate  is 
clamped  to  a  rigid  foundation  (the  plane  z  =  (  —  h) 


Because  of  symmetry  only  1  /4'th  of  the  domain  is  analyzed.  Material  properties 
are  assumed  isotropic  and  linearly  elastic.  The  modulus  of  elasticity  is  denoted  E 
and  Poisson'.s  ratio  Cartesian  displacements  and  stresses  are  denoted  u,  n,  w  and 
Gj.Oy.a^.  Tj-y.  Tj,,  and  Ty,  respectively. 
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The  face  ABC D  is  loaded  by  the  parabolic  tangential  shear, 

Trz  =  4(z/0(l-z/0 

Txy  =  0 


where  t  is  the  thickness. 

The  T-dispiacements  are  u=0  on  the  face  ABCD.  The  face  IJKL  is  clamped 
i.e. 

u  =  V  =  U’  =  0  (2) 

On  the  face  BFKJGC  the  symmetry  conditions  are, 

r  =  0 

V  =  0  (3) 

Txy  =  0 

All  other  faces  are  traction  free. 

This  problem  was  analyzed  as  a  benchmark  problem  in  [11]  and  [12].  It  is 
an  excellent  academic  benchmark  problem  illuminating  the  difficulties  of  reliable 
computations.  We  will  address  here  the  dependence  of  cri.(T,0,0)  on  Poissons  ratio 
1/  and  the  behaviour  of  the  solution  close  to  the  vertex  E  and  the  edge  EF. 

Let  us  list  the  major  theoretical  properties  of  the  solution  to  this  mathematical 
problem. 


•  There  exists  exactly  one  solution  which  has  finite  energy 

•  The  solution  is  analytic  everywhere  on  Q  and  its  faces  except  for  the  edges. 
Because  of  symmetry  the  solution  is  also  analytic  on  the  edge  JGCBFK 

•  The  solution  is  singular  along  edges  and  the  vertices 


In  a  neighbourhood  of  the  edges  and  the  vertices  the  solution  can  be  written  as  a 
sum  of  singular  terms  and  a  smoother  function.  The  singular  behaviour  is  different. 

•  along  the  open  edges 

•  in  a  neighbourhood  of  the  vertices 


2.1  Solution  behaviour  close  to  edges 

We  will  discuss  the  properties  of  the  exact  solution  close  to  the  edge  EF. 
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Denoting  by  r  tb  distance  of  the  point  (x,y,z)  from  edge  EF  and  by  6  (0  < 
0  <  37r/2)  the  pola’  ;gle,  the  displacements  have,  for  any  s  >  1  and  s 
the  form  [13] 


n(r,  6.y) 

.  U'(r,  e*,!/) 

+  E  (4) 

+  E  A'j;f +  smoother  term 


flel-X /c<3 


where  for  A'  zero  or  even  (a  =  /,  II  or  III  in  the  following)  are  of  the  form 


■  ■ 

r 

-5- 

II 

0 

II 

1 

and  for  k  odd 

0 

II 

0 

0 

0 

^n.3  (^)  . 

0 

0 

,  ^\f,\e]  = 

■  ■ 

0 

0 

^1,13(0) 

The  edge  intensity  functions  A’[,''*(i/)  for  A-  >  0  can  be  expressed  as  a  linear 
combii,ation  of  derivatives  of  K[;'^\y  ].  In  the  sequel  we  will  use  the  notation  K''J'[y) 
=  E]J'^'\y].  In  technical  applications  A/}’()y)  and  K]]](y)  are  the  import atU 

edge  intensity  functions. 

Here  functions  depend  only  on  local  geometry,  material  and  local  boundar}' 

conditions.  In  fact  the  functions 


r^i" 

■  ■ 

0 

A?.’ 

,  r 

■  ■ 

0 

■  1 

3 

are  identical  with  the  singular  functions  for  two-dimensional  elasticity  and  the 
functions 

0 


r  III 


’1.2 

0 
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are  the  singular  functions  for  the  torsion  problem,  i.e.  the  Laplace  operator. 

Edge  stress  intensitt  functions  K\]\y)  and  K\]\{y)  are  analytic  on  the 

semi-open  interval  0  <  y  <  u’/2,  singular  at  y=0  and  depend  on  the  global  data  of 
the  problem. 

The  values  A*,-"’  can  also  be  complex.  In  this  case  the  functions  are 

complex  too  and  the  real  part  of  (4)  is  taken  instead.  In  this  case  the  solution  has 
an  oscillatory  character.  In  general,  for  exceptional  cases  f geometrical  angles  or 
material  data),  besides  the  term  there  can  also  be  terms  of  the  form  r''(logr)’ 
(we  will  not  address  this  Cctse  here). 

In  (4)  we  assume  that  the  singularity  exponents  are  numbered  such  that  < 

/?e[ A),’ ■*■’*]  .  The  stresses  are  infinite  if  any  Ae[AW]  <  1  and  tne  corresponding  edgt- 
intensity  function  is  nonzero. 


A</' 

Aj  .  A, 

x(2)  W3) 

A,j  .  All 

0.544 

1.629  ±0.231? 

0.909 

2.301  ±  0.316? 

Table  1:  First  edge  eigenvalues  A/’*  and  A*//  for  edge  EF 


a‘" 

a'*' 

\(3) 

^ni 

0.GG7 

1 

1.333 

2 

Table  2:  First  edge  eigenvalues  A;;/  for  edge  EF 


I'or  the  model  example  the  first  edge  eigenvalues  A*/'  are  given  in  tabh^s  1  and 
2.  These  values,  which  arc  independent  of  Pois.'ion's  ratio  arc  obtained  by  solvim; 
a  nonlinear  equatioii  [l-l|.[l5p  We  see  that  the  first  ternt  for  the  mode  /  .  mode  I  / 
aiid  mode  III  components  of  (4)  yield  infinite  stresses  at  the  edge  EF. 

1  he  edge  stress  intensity  function  A’|''*(y)  has  in  a  neighbourhood  of  y  =  0  the 
following  form 


^  •''Wft"''’ ‘  +  smoother  term  (T 

/=! 

The  intensity  factor  is  the  mode  I  vertex-edge  intensity  factor  of  order  j  for 
edge  EE  at  vertex  E.  The  exponent  is  related  to  the  values  of  A*/*  and  .A*/* 
characterizing  the  edge  and  the  vertex  singularities  as  discussed  below.  .Analogous 
expressions  apply  for  the  mode  11  and  mode  JIl  edge  stress  intensity  functions. 
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2.2  Solution  behaviour  close  to  vertices 


In  a  neighbourhood  of  the  vertices  E,  L  and  /  the  solution  ha.s  a  singular  behaviour 
which  in  spherical  coordinates  can  be  expressed  as 


u(r,^\0) 

j 

■  0l^  Vs  0)  ■ 

v{r,u.\  0) 

-f  smoother  terms 

tr(r.  w^-,  0) 

j=i 

(6) 


In  (6)  the  coefficients  5*-''  are  the  rcr.ez  intensity  factors.  The  functions  0, 
caii  be  understood  to  be  defined  on  the  intersection  of  the  domain  and  a  spherical 
surface.  This  intersection  has  a  boundary  with  the  corners  (LJk,<Pk)  associated  to 
the  edge  k  of  the  domain. 

For  angle  (^'k.  Ok )  associated  with  edge  EF a.nd  for  i/  >  0  the  functions  0,*^'.  0 
(with  ?  =  1,2,, 3)  are  found  (in  the  numerical  analysis)  to  correspond  to  a  pure  mode 
1  deformation  at  the  edge  EF.  The  functions  0,  and  0,  '  correspond  to  a  combined 
mode  II  and  mode  III  deformation  of  EE  with  a  zero  mode  7  component. 


1/ 

,\(i) 

A(3) 

A(4-5) 

0.0 

0.544 

0.909 

1.054 

1.629 

0.30 

0,025 

0.785 

1.237 

1.521  ±  0.172! 

0.43 

0.711 

0.751 

1.310 

1.331 

iable  3:  The  \ertcx  coefficients  .\F)  for  vertex  E 


For  vertex  E  the  first  vortex  singularity  exponents  are  given  in  table  3.  By 
(■o::i;)arii;g  (4i  and  (61  we  get  the  coefficients  in  (5)  as 

=  .v"  -  A'/'.  (7) 

In  (7,i  we  consider  only  /-values  which  are  associated  with  a  mode  I  deformation 
of  edsc  EF  {( .g.  I  =  I  5. . . 


1  =  1 

/  =  2 

/  =  3 

/  =  4  and  /  =  5 

9i.t 

0.081 

- 

- 

0.977  ±  0.172! 

-.0) 

^II.i 

- 

-0.123 

0.329 

- 

'III.i 

- 

0.119 

0.571 

- 

Table  ■}:  The  coefficients  and  i  for  edge  EF  at  vertex  E  for  the  case 

i'=0.3 


In  table  4  the  first  coefficients  are  given  for  the  case  1/  =  0.3. 
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V\  e  see  that  edge  stress  intensity  functions  can  be  bounded,  go  to  zero  or  can  be 
unbounded  for  y  small.  We  mention  that  in  our  example  some  in  (5)  could  be 
zero.  In  fact  this  occurs  in  the  case  i/  —  0  because  the  solution  is  t/-independent. 

On  the  edge  AE  in  a  neighbourhood  of  the  vertex  E  we  get,  using  (6),  the  stress 

j 

Ox  =  —  +  smoother  term.  (8) 

1=1 

The  stress  Ox  is  singular  at  vertex  E  for  all  values  u  provided  that  the  coefficients 
9i  0  or  92  ^  0.  As  we  will  see  later  for  particular  values  of  u,  some  coefficients  9, 
could  be  zero  (especially  91)  and  hence  it  is  essential  that  two  terms  in  (6)  and  (S) 
are  determined. 

We  have  briefly  described  the  properties  of  the  exact  problem  to  our  benchmark 
problem.  For  a  theoretical  discussion  of  the  convergence  properties  of  the  p-version 
in  three  dimensions,  see  ([16], [17]). 

The  characterization  of  the  three-dimensional  solution  is  of  great  importance 
since  the  properties  of  the  exact  solution  influence  drastically  the  performance  of 
the  finite  element  method  and  have  direct  impact  on  the  proper  design  of  the  mesh,. 
Let  us  note  that  for  the  design  of  a  mesh  we  do  not  have  to  know  the  exact  character 
of  the  singularities.  In  fact  we  normally  do  not  know  them.  We  nevertheless  know 
the  qualitati\e  character  of  the.se  singularities  and  we  can  design  a  good  mesh  fol¬ 
lowing  simple  a-priori  known  rules.  The  mesh  which,  if  optimal,  should  be  refined 
(cylindrically)  in  a  neighbourhood  of  edges  and  in  a  neighbourhood  of  the  vertices 
(spherically). 

Let  us  now  discuss  some  numerical  results  for  the  following  model  dimensions. 

b/t  =  300 

h/i  ^  50  (9) 

u-/t  =  100 

t  =  100 

w  e  will  present  results  computed  by  the  finite  element  program  STRIPE  based 

oil  a  />-\ersion  of  tlie  finite  element  method.  Results  obtained  with  the  /i-version  of 
FE.M  are  given  in  [l  1  ]. 

In  figure  2  we  show  the  mesh  used  for  the  computation.  In  a  neighbourhood  of 
the  edge  EE  and  the  vertex  E  we  used  a  strongly  refined  mesh  as  shown.  The  mesh 
used  ha.s  three  spherical  layers  of  elements  around  vertex  E  and  six  cylindrical  layers 
of  elements  around  edge  EE.  Using  this  mesh  in  combination  with  the  p-version  of 
FE.M,  reliable  results  will  be  obtained.  The  mesh  used  has  262  elements. 

In  table  5  we  report  the  energy  11^111(0)  of  the  finite  element  solution  as  a  function 
of  p  and  1/  and  the  relative  error  ||c||£’(n)  =  ||u  -  u||£(n)  the  energy  (the  exact 
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p 

1/  = 

=  0 

ll'^l(£(n) 

1/  — 

0.30 

Ikilltm 

0.45 

lkil£:(m 

2 

5.60000S 

3.500017 

4.959967 

3.680763 

3.480780 

4.584843 

3 

9.0540S9 

0.045936 

8.414659 

0.226071 

7.563369 

0.502254 

4 

9.09SS95 

0.001130 

8.600359 

0.040371 

7.961862 

0.103761 

5 

9.099S61 

0.000164 

8.638234 

0.002496 

8.058791 

0.006832 

6 

9.099963 

0.000062 

8.639845 

0.000885 

8.063074 

0.002549 

7 

9.099999 

0.000026 

8.640380 

0.000350 

8.064567 

0.001056 

Estim. 

9.100025 

8.640730 

8.065623 

Table  5:  The  energy  ||i/||£(n)  energy  error  of  the  FE-solution  ||e|||:(Q)  ob¬ 

tained  with  the  262  element  mesh.  Data  given  are  scaled  with  a  factor  10“’^ 

energy  was  estimated  by  extrapolation).  We  see  that  the  p-version  provides  a  very 
good  assessment  of  the  accuracy  achieved  and  that  the  accuracy  is  influenced  by  the 
Poisson  ratio 

Let  us  note  that  the  convergence  with  respect  to  p  has  two  phases.  The  first  one 
(preasN’mptotic)  has  exponential  convergence  rate  while  the  second  has  algebraic 
convergence  rate.  The  reason  is  that  for  high  p  the  mesh  is  ’’underrefined"  and  that 
the  rate  is  algebraic  (see  [IS]  and  [19])  because  the  governing  error  resides  in  the 
elements  which  contains  the  edges  and  vertices.  The  transition  between  the  two 
phases  occurs  roughly  when  errors  in  these  elements  becomes  the  dominating  part. 
The  design  of  the  mesh  is  near  optimal  when  the  required  accuracy  is  achieved  at 
the  end  of  the  transition  part. 

In  general  some  overrefinement  is  recommended  because  we  would  like  to  avoid 
having  an  algebraic  rate  of  convergence.  The  mesh  in  a  neighbourhood  of  the  sin¬ 
gularity  is  recommended  to  be  geometric  with  the  ratio  0.15.  This  ratio  is  in  some 
sence  optimal.  For  a  theoretical  analysis  of  the  one-dimensional  model  which  is 
rouglily  ap[)licable  here,  see  [20],[2l|. 

Figure  3  shows  in  a  lin-log  scale  the  stress  0, 0)  and  C7i(x.  u'/2, 0)  for  various 
values  of  the  Poisson  ratio  i/.  Stresses  have  been  obtained  by  direct  calculation  from 
the  displacements.  For  x  very  small  there  is  a  considerable  jump  in  computed 
stresses  at  element  boundaries.  We  see  that  at  the  edge  AE  the  stress  distribution 
is  ver}  different  from  the  stress  at  the  line  (x,i(’/2,0).  For  n  =  0  there  is  no  difference 
between  these  behaviours.  Close  to  vertex  E,  we  see  the  strong  influence  of  Poisson's 
ratio  on  the  stress  . 

Table  6  to  Table  8  shows  the  coefficients  and  q2  and  vertex  exponents  A'''  and 
computed  ([8], [22])  using  the  mesh  shown  in  figur  b  The  coefficients  q,  have 
been  obtained  from  computed  A-values  and  eigenfunctions  (see  (6)). 

Values  labelled  ’’exact"  in  table  6  have  been  obtained  with  a  very  detailed  two- 

dimensional  model  and  p=10  uniformly.  All  digits  of  the  "exact”  values  are  believed 


16 


FFA  TN  1992-17 


Figure  3:  Calculated  stress  distributions  (7i(x,0,0)  and  (7i(x,  u-72, 0)  for  benchmark 
problem  (p  =  7  solution) 
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p 

Ad) 

A(2) 

9i 

92 

3 

0.5442 

0.9079 

-73.06 

14.67 

4 

0.5444 

0.9080 

-75.24 

6.11 

5 

0.5445 

0.9085 

-75.81 

-1.39 

6 

0.5445 

0.9085 

-75.87 

-1.55 

7 

0.5445 

0.9085 

-75.87 

-1.61 

Exact 

0.5445 

0.9085 

-75.89 

-1.61 

Table  6:  Calculated  A-values  and  ^-values  as  function  of  p  for  i/  =  0 


P 

Ad) 

\{2) 

9i 

92 

3 

0.62.56 

0.7871 

-10.74 

-1.10 

4 

0.6254 

0.7853 

0.06 

-9.77 

5 

0.62.53 

0.7851 

-1.13 

-9.01 

6 

0.62.53 

0.7852 

-0.45 

-8.70 

7 

0.62.53 

0.78.52 

-0.22 

-8.68 

Table  7:  Calculated  A-values  and  g-values  as  function  of  p  ior  u  =  0.30 


P 

Ad) 

A(2) 

9i 

92 

3 

0.7115 

0.7548 

26.51 

-2.70 

4 

0.7118 

0.7514 

44.43 

-11.90 

5 

0.7116 

0.7510 

42.45 

-10.35 

6 

0.7116 

0.7511 

43.54 

-9.84 

7 

0.7116 

0.7511 

43.80 

-9.84 
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tc  be  significant. 

W'e  see  that 

a)  the  accuracy  of  computed  data  can  be  judged  from  the  sequence  of  solutions 
generated 

b)  the  coefficients  A  and  q  can  be  accurately  computed  (except  for  qi  when  i/  = 
0.3) 

c)  for  =  0.3  the  first  vertex  intensity  factor  5^*^  occasionally  is  very  small 
leading  to  a  very  small  region  with  high  stresses  (compare  figure  3) 

d)  for  2'  >  0.3  we  have  the  paradoxical  result  that  for  x  small  the  stress  is 
infinite  with  sign  opposite  to  what  one  would  expect 

Figure  4  show  the  energy  error  ||e  lictn)’  the  error  in  values  and  Qi  for 

the  case  r  =  w  e  see  that  the  error  in  A*’*  and  in  q\  has  same  character  as  the 
accuracy'  of  the  energy. 

In  figure  5  we  show  in  lin-log  scale  the  graph  of  the  function  (7j.(  0,  0)  as  function 

of  I  for  various  Poisson's  ratios  when  computed  from  the  asymptotic  expansion  (8) 
using  the  first  7  terms.  The  asymptotic  expansions  are  in  good  agreement  with 
stre.'^.'^es  calculated  directly  from  the  FE-solution  with  p  =  7. 

Let  us  now  address  the  stresses  at  the  edge  EF.  Here  stresses  are  infinite  and  we 
have  to  compute  the  edge  stress  intensity  functions  K\]\y)  and  K 

in  {y)-  In 

a  forthcoming  paper  [10]  we  will  describe  the  procedure  used.  Briefly,  the  intensity 
functions  are  approximated  by  piecewise  polynomials  (or  trigonometric  functions) 
a.' 

~  Xc,Q,(y)  (10) 

1=0 

where  Q:[y)  is  the  Legendre  polynomial  of  order  i  and  c,  are  coefficients  to  be 
calculated.  W'e  compute  the  coefficients  c,  when  EF  is  divided  into  the  three  sub- 
interval-s  i)  <  y  <  t A  <  y  <  It.  It  <  y  <  50<  (compare  figure  2). 

In  table  9  we  show  the  calculated  values  of  A^j'*(j/)  (from  (10))  in  the  three 
mentioned  intervals  for  2/  =  0.3  as  function  of  p.  Further,  A'j’\0)  is  found  to  be 
close  to  zero  as  predicted  from  (5)  and  table  4.  For  higher  values  of  p  the  calculated 
intensity  function  A’j"(p)  is  continuous  over  the  interval  0  <  y  <  w/2. 

In  figure  6  we  show  in  log-log  scale  the  function  K\^\y)  (the  p  =  8-solution)  for 
1^  =  0. 3  and  x  small.  Since  qi  is  almost  zero  (table  7)  we  expect  that  the  second  term 
in  {')}  influencing  A’}'’(y)  might  be  the  dominating  one  (except  for  y  very  small). 
Figure  G  reveals  that  this  is  the  case.  Note  that  because  and  are  complex 
A'l’^(y)  is  oscillatory  but  only  in  a  close  neighbourhood  of  y  =  0. 
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Figuro  4;  F2rror  in  calculated  data 
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Domain 

y' 

p  =  8 

p  =  l 

p  =  6 

p  =  5 

p  =  4 

p  =  3 

p  =  2 

1 

0 

-1 

-1 

-6 

-24 

-9 

-225 

-421 

1 

tjb 

-158 

-1.59 

-161 

-173 

-151 

-407 

-628 

1 

2t/.5 

-328 

-329 

-331 

-344 

-320 

-574 

-797 

1 

'it  lo 

-1S2 

-482 

-482 

-493 

-468 

-720 

-930 

1 

4t/5 

-615 

-615 

-614 

-622 

-585 

-840 

-1026 

1 

t 

-730 

-728 

-724 

-730 

-698 

-931 

-1086 

2 

0 

-730 

-729 

-723 

-732 

-691 

-936 

-1114 

2 

fd/.o 

-1146 

-1145 

-1143 

-1147 

-1148 

-1402 

-1432 

2 

12t/5 

-1331 

-1328 

-1322 

-1320 

-1323 

-1598 

-1648 

2 

lSl/5 

-142S 

-1424 

-1415 

-1403 

-1370 

-1622 

-1761 

2 

24t/5 

-1485 

-1481 

-1469 

-1444 

-1384 

-1570 

-1772 

2 

6/ 

-1522 

-1516 

-1499 

-1463 

-1405 

-1541 

-1681 

3 

0 

-1522 

-1516 

-1499 

-1461 

-1405 

-1533 

-1675 

3 

Aitjb 

-15S8 

-1591 

-1.599 

-1620 

-1627 

-1571 

-1603 

3 

S6t /b 

- 1 58 1 

-1578 

-1576 

-15.58 

-1593 

-15.54 

-1554 

3 

2^M;b 

-1569 

-1571 

-1576 

-1566 

- 1 535 

-1508 

-1527 

3 

\12t/b 

- 1 55S 

-1.561 

-1554 

-1582 

-1563 

-1458 

-1523 

3 

43/ 

-1.561 

-1547 

-1.566 

-1531 

-1664 

-1429 

-1541 

lalilf  9:  Calculated  values  of  the  mode  /  edge  intensity  function  at  six 

equidistant  points  inside  domains  1  -  3  for  =  0.3.  Plate  thickness  is  (  and  local 
y-coordinatc  inside  subdomains  is  y' 
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In  figure  7  we  depict  the  pointwise  accuracy  of  A'j’*(y)  for  0  <  y  <  t  and  the 
accuracy  of  the  energ}'.  As  exact  values,  the  solution  for  p=S  is  used.  Once  more  v.  e 
find  that  the  accuracy  of  the  computed  pointwise  values  of  Kj^\y)  are  essenliallv 
the  same  as  the  accuracy  of  the  global  energy. 

2.3  The  stress  concentration  problem 

In  practice  the  designer  tries  to  avoid  sharp  (internal)  edges  and  vertices  as  to 
prevent  fatigue  crack  initiation.  Relevant  design  parameters  for  nucleation  of  such 
cracks  are  the  magnitude  of  the  local  stress  components  and  their  time  history  [2\ 
The  designer  has  to  have  an  effective  tool  to  get  the  relationship  between  stress  and 
radius  at  a  fillet.  By  using  asymptotic  analysis  (cis  in  the  case  of  the  analysis  of 
stresses  in  a  neighbourhood  of  sharp  edges  and  corners)  a  rational  design  formula 
can  be  developed. 

Consider  the  solution  in  a  neighbourhood  of  the  radius  Ry  (figure  8)  and  assume 
for  a  moment  (for  reasons  of  clearity)  that  the  solution  is  two-dimensional.  For  a 
sharp  edge,  i.e.  R\  =  0,  the  solution  may,  for  any  s  >  1  and  s  /  *]■  be  written 

u{r.O}  ' 

v{r,6)  =  ^  -f  smoother  term  (11) 

_  u'ir.O]  \  Re[\‘J']<s 

where  are  the  singular  functions  of  two-dimensional  elasticity  and  the 

Laplace  equation.  A'),-^*  are  constants. 

We  arc  interested  in  solutions  for  i?i  >  0  in  a  small  region  (of  size  2R\.  say> 
around  the  (rounded)  edge  EF.  W'e  first  give  the  asymptotic  behaviour  of  the 
displacements  for  R\  tending  to  0. 

These  displacements  are  characterized  by  the  functions  ^1^3(1-,;)  defined  on  an 
infinite  two-dimensional  domain  with  a  rounded  vertex  with  radius  R  =  \  (figure  9 
with  r  — »  c>c  and  w'  =  r/2  ).  Denote  by  the  stresses  corresponding  to  the 

displacements  nj]^j(x.r).  The  functions  u|'’^(x,r)  are  solutions  to  the  equations  of 
two-dimensional  elasticity  with  boundary  conditions  that  displacements  are  of  the 
type  ^[/>(0)  for  r  large. 

The  stresses  at  a  point  (x.s)  in  a  neighbourhood  of  the  radius  Rj  may  then  be 
written 


(7fc((x,z)  = 


(12) 


T 
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wliere  x  =  x/Ri  and  5  =  :IR\- 

Here  =  min_,,Q  /?e[AW]  is  the  singularity  exponent  with  the  smallest  real  part, 
and  (12)  gives  the  asymptotic  behaviour  for  (i,i)  fixed  and  Ri  going  to  0  [23]. 
Therefore  we  can  expect  that  the  sum  in  (12)  gives  a  good  approximation  for  the 
stress  near  the  rounded  edge  if  is  sufficiently  small  (i.e.  ’■elative  to  the  plate 
thickness). 

For  the  three-dimensional  case,  the  full  mathematical  theory  is  not  available. 
Nevertheless,  for  R\  going  to  0,  we  have  the  following  expression  for  the  stresses  at 
a  point  {x.y.z)  close  to  edge  EF 


(7k:{x.y.:) 


H  ’  +  higher  order  terms 

a=i.n.jii 


(13) 


In  practice,  the  functions  are  determined,  once  for  all,  by  solving  two 

ela.sticity  problems  and  the  Laplace  equation  on  a  two-dimensional  domain  (figure  9). 
For  three-dimensional  problems,  we  only  have  to  compute  the  solution  on  the  domain 
with  sharp  edges  and  extract  the  vertex  intensity  functions  A'^’*(i/).  Then  we  can 
immediately  predict  the  stresses  for  small  Ri  >  0  using  (13). 

We  here  give  numerical  results  from  an  analysis  of  the  folded  plate.  The  case 
r  =  0  is  considered.  For  ^  =  z  j  2  we  are  interested  in  A’j'*(t/),  K\'\y).  K\^\y)  and 
A’]|  (;;)  while  the  corresponding  edge  eigenvalues  satisfy  <  3p  with  p  —  X'/'. 


p 

Af> 

A’)^*  and 

Rr 

A'd) 

2 

-150.5.9 

-0.525 

± 

1.278/ 

-7.85 

3 

-145T3 

-0.3G8 

± 

2.041/ 

-7.19 

4 

■1482.9 

-2.445 

X 

5. 4  or./ 

-7.90 

5 

-1490.2 

-2.495 

± 

5.450/ 

-7.93 

ti 

•1490.8 

-2.511 

± 

5.514/ 

-7.91 

( 

-1491.0 

-2.513 

± 

5.520/ 

-7.91 

8 

-1491.1 

-2.514 

± 

5.524/ 

-7.91 

9 

-1491.1 

-2.514 

± 

5.-525/ 

-7.91 

10 

-1491.1 

-2.514 

± 

5.525/ 

-7.91 

Table  10:  Calculated  edge  functions  A'^-'^  for  u  =  0 

Table  10  shows  calculated  edge  intensity  functions  obtained  with  a  mesh  having 
six  cylindric  al  layers  of  elements  (grading  factor  0.15)  around  edge  EF  (radius  Ri  = 
0).  Table  11  sliows  the  relative  error  in  estimated  maximum  von  Mises  effective 
stress  for  tlirc-e  different  Aj -values.  The  error  is  defined  as  the  estimated  value  (from 
ecjuation  ( 12j )  minus  the  exact  value.  As  exact  solutions  we  have  used  FE-solutions 
obtained  with  p  =  12  uniformly  on  very  detailed  two-dimensional  meshes. 
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p 

II 

o 

o 

Rx/t  =  0.1 

Rx/t  =  0.25 

2 

0.015 

0.06 

0.11 

3 

0.020 

0.03 

0.10 

4 

0.002 

0.05 

0.17 

5 

0.003 

0.06 

0.17 

6 

0.004 

0.06 

0.18 

7 

0.004 

0.06 

0.18 

8 

0.004 

0.06 

0.18 

9 

0.004 

0.06 

0.18 

10 

0.004 

0.06 

0.18 

Table  11:  Maximum  relative  error  in  von  Mises  effective  stress  at  fillet.  Maximum 
effective  stresses  are  2224,  753,  and  480  units  for  i?i/t=  0.01.  0.1  and  0.25  units 
respectively 

The  results  in  table  11  shows  that  the  estimate  (12)  is  very  accurate  for  Ri 
small.  The  main  error  of  the  method  in  this  case  is  (already  for  p  >  4)  caused  bv 
the  term  and  not  the  discretization  error.  From  equation  (12)  we  see  that 

the  rtlatu'i  error  in  stress  behaves  like  for  Rx  tending  to  0.  With  p  =  0.544 
(table  1 )  we  expect  that  the  relative  error  should  behave  like  Data  in  table  1 1 

show  that  this  is  the  case. 

For  the  benchmark  problem  considered,  we  may  simplify  (12)  somewhat  since 
A/"  by  far  is  the  largest  intensity  factor  (table  10).  By  truncating  the  series  (12) 
to  one  term  the  following  expression  for  the  maximum  von  Mises  stress  at  the  fillet 
is  obtained. 

(14) 

where  C’  =  275. 

Ihis  expre.^sion  is  asymptotically  exact  for  A  — ►  0  (with  tree  digits  accuracy) 
while  for  Rjt  —  0.01.  Rjt  =  0.10  and  Rjt  =  0.25  the  relative  errors  in  maximum 
stress  are  0.77(.  4.277  and  7.7%  respectively. 

A  similar  approach  as  shown  here  for  edges  can  be  used  for  analysis  of  rounded 
vertices. 


I  he  procedure  developed  makes  it  possible  to  efficiently  and  reliably  compute 
poiutwif'C  stresses  while  the  error  in  edge  intensity  functions  converge  as  the 
error  in  energy.  Since  (13)  is  an  ana/y/tca/ expression  for  the  stresses  at  the  fillet,  it 
is  suital>le  for  use  in  complex  dimensioning  formulas  (fatigue  crack  initiation  etc.). 
W  e  note  that  only  ont  finite  element  analysis  is  needed  (with  7?  =  0)  in  order  to 
derive  the  unknown  functions  A'^-’*  in  (13). 
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3  THE  FATIGUE  DESIGN  OF  COMPLEX 
THREE-DIMENSIONAL  COMPONENTS 

In  this  section  we  demonstrate  the  reliable  analysis  of  stresses  at  fillets  in  a  complex 
aircraft  fuselage  frame  (figure  10)  using  the  techniques  reviewed  in  section  2. 

We  analyze  models  having  0.1  -  1.2  million  degrees  of  freedom.  Our  global  finite 
element  model  (designed  for  the  /i-version  of  the  finite  element  method)  hcis  5-510 
elements.  The  sub-section  to  be  used  to  demonstrate  the  reliable  analysis  of  stresses 
and  stress  intensity  factors  using  h—p  type  meshes  has  600-800  elements  (depending 
on  the  goals  of  the  analyses). 

We  assume  that  the  material  is  linear  isotropic  and  homogenuous  with  Poisson's 
ratio  i/=0.3.  Traction  boundary  conditions  applied  on  the  sub-section  considered 
are  obtained  from  an  analysis  of  the  global  model  (figure  10a). 

3.1  Stresses  at  hole  boundaries 

The  domain,  as  in  the  section  2,  has  several  smooth  closed  edges  (labelled  A.B....). 
They  are  similar  to  the  edge  .AE  in  figure  1.  The  main  difference  is  that  edges  are 
curved.  The  cur%ature  of  the  edges  complicates  the  form  (4)  for  i  >  1.  Additional 
terms  of  the  type  r'^'*  ■*’*''I'^-’'*'*(0)(with  k  integer  and  ^  >  1)  have  to  be  added  to 

(4j.  For  edges  A -  we  have  =  2.740  -f  1.119r,  =  4. SOS  A  1.464?  and 

'^/V/  =  Because  /f€(Ao)  >  1,  the  stresses  are  essentially  smooth  and  hence  the 
stress  intensity  functions  are  not  of  interest. 

Because  of  assumed  boundary  conditions  (free  of  tractions  in  a  neighbourhood 
of  the  holes),  the  only  nonzero  stress  at  the  edge  is  the  normal  stress  along  the  edge. 
If  the  shape  of  the  edges  is  analytic  then  the  stresses  are  analytic  functions  in  the 
edge  length  parameter.  This  function  is  of  special  interest  since  often  the  stresses 
have  local  maxima  at  the  hole  boundaries. 

We  essentially  proceed  as  in  section  2  when  we  computed,  in  a  postprocessing 
mode,  the  coefficients  of  the  expansions  in  the  Legendre  polynomials  (10).  In  a 
similar  way  we  can  compute  stress  components  a,j  on  smooth  surfaces  by  using 
the  double  sum 


p  p 

2/2)  ~ 


1=0  }=0 


(15) 


where  yi  and  j/2  are  surface  coordinates.  The  advantage  of  computing  the  coef¬ 
ficients  c,j  and  using  (15)  for  stress  calculation  is  that  the  errors  in  coefficients  c,j 
decrease  as  the  energy  error.  The  problem  of  computation  of  stresses  along  smooth 
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h-version  mesh 


Figure  10:  Fuselage  frame  (a)  and  sub-section  for  detailed  analysis  (b) 
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edges  and  on  smooth  surfaces  we  address  in  [23]. 


3.2  Stresses  at  fillets 

Below  we  elaborate  on  some  practical  aspects  of  computations  of  the  problem  of  the 
fuselage  shown  in  figure  10. 

The  mesh  has  to  be  designed  so  that  the  error  of  the  Fli-solution  decrease  ex¬ 
ponentially  with  increasing  number  of  degrees  of  freedom  N  in  the  p— range  used. 
We  have  at  our  disposition  (section  2)  postprocessing  techniques  which  allows  us  to 
compute  engineering  parameters  of  interest,  that  is  stresses,  stress  intensity  factors 
etc.  with  the  rate  of  convergence  as  the  energy’  (not  energy  norm). 

When  employing  the  h  —  p  version  of  the  finite  element  method,  the  optimal 
mesh  should  be  very  strongly  graded  in  a  proper  sense  towards  the  edges  and  ver¬ 
tices  (compare  figure  2).  With  known  characteristics  of  the  solution  (eg  location  of 
boundary  layers,  strengths  of  singularities  etc.),  it  is  possible  using  few  simple  rules 
for  mesh  design  to  obtain  high  accuracies  and  an  effective  control  of  the  error.  A 
mesh  generator  Wcls  implemented  in  the  FE-code  STRIPE  which  uses  a  coarse  mesh 
(prepared  by  the  user)  to  create  a  mesh  strongly  refined  in  a  neighbourhood  of  edges 
and  vertices. 

Figure  11a  shows  the  coarse  user  defined  FE-mesh  {R,  =  0,  t  =  1,2, . .  .4)  mod¬ 
elling  the  sub-section  of  the  fuselage  (figure  10b).  The  generated  mesh  for  the  h  —  p 
version  is  shown  in  figure  11b.  We  have  generated  additional  element  layers  at 
internal  edges  only  (where  <  1). 

Blended  function  mapping  is  used  [24]  to  exactly  model  the  shape  of  the  domain 
(only  the  polyhedral  shape  is  depicted  in  figure  11).  Note  that  the  details  of  the 
generated  mesh  are  invisible  at  this  scale.  The  input  data  for  the  mesh  generator  is. 
except  the  coarse  mesh,  a  list  of  the  edges  and  vertices  where  the  mesh  should  be 
refined  and  the  number  of  layers  of  elements  to  be  generated  around  the  vertices  and 
edges.  The  h  —  p  approach  taken  has  the  advantage  that  the  global  mesh  generation 
problem  splits  into  a  number  of  local  mesh  refinement  problems.  A  mesh  generator 
based  on  these  principles  is  relatively  simple  to  implement. 

Figure  12  shows  the  calculated  von  Mises  effective  stress  distribution  in  the  frame 
for  p=6.  The  stress  distribution  shown  is  very  close  (with  actual  resolution)  to  the 
exact  solution.  The  figure  shows  that  for  actual  radices  R,  =  IfSmm  analyzed  high 
local  stresses  at  the  edges  are  obtained.  The  edge  labelled  Ri  which  seems  to  be 
most  critical  will  be  studied  in  detail  below. 

Let  us  address  the  problem  of  determining  the  smallest  radices  R  (i  =  1, 2, ...  4) 
such  that  the  maximum  stress  at  the  stress  concentration  R\  ,  for  example,  will  be 
the  same  a.s  the  one  at  the  edge  B  (figure  10b).  This  is  done  by  the  technique 
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described  in  section  2.3.  First  we  compute  the  stress  intensity  functions 
K\]\y)  and  along  the  unrounded  edge  and  from  (13)  we  compute  the 

stresses. 

Calculated  edge  stress  intensity  functions  for  different  degrees  p  of  the  el¬ 
ements  along  the  edge  are  shown  in  figure  13.  Three  layers  of  elements  have  been 
generated  around  each  edge  and  vertex  (figure  11).  Three  subdomains  along  the 
edge  are  used  for  extraction  of  the  polynomial  coefficients  c,  (10).  The  results  in 
figure  13  show  for  p=5  a  smooth  and  continuous  variation  along  the  edge.  With 
increasing  p  the  edge  intensity  functions  quickly  converge. 

Our  experience  is  that  by  using  four  layers  of  elements  with  the  grading  factor 
0.15.  an  exponential  rate  of  the  convergence  is  obtained  in  the  range  3  <  p  <  7 
with  approximat  .ly  the  error  ~  C  ■  10“°®'’  (for  energy,  stress  intensity  factors  etc.) 
provided  that  the  mesh  in  the  entire  domain  is  reasonable.  In  3D  situations  the 
global  mesh  often  becomes  partly  distorted  or  element  aspect  ratios  become  ver}' 
large  (for  simplicity  of  input  data  generation)  which  may  slow  down  the  convergence 
rate  for  lower  p-levels. 

Figure  14  shows  the  maximum  von  Mises  stress  at  any  point  ((j.  c)-section)  in  the 
fillet  as  function  of  the  position  y  along  the  edge.  Stresses  have  been  calculated  using 
(13)  and  stress  intensity  data  shown  in  figure  13.  For  comparison,  stresses  obtained 
from  a  detailed  FE-analysis  with  p  =  6  (figure  12)  with  radius  R\  =  lj3mm  are 
shown.  Excellent  agreement  between  estimated  (13)  and  calculated  local  stresses 
is  obtained.  Stresses  estimated  from  (13)  for  different  /?i-values  are  also  shown  in 
figure  14.  With  radius  /?i  =  3.0mm  the  maximum  stress  in  the  fillet  becomes  tlie 
same  as  the  maximum  stress  at  the  hole  boundary  (460  Mpa.  see  figure  12). 

To  conclude,  the  use  of  mesh  generation  in  a-priori  known  areas  where  solution 
exhibit  strong  boundary  layers,  use  of  postprocessing  techniques  to  calculate  A’0>(y ) 
in  (13)  gives  a  rational  design  tool  for  sizing  of  fillets  in  complex  three-dimension;!! 
components. 


Kl(y).  Kij(y),  Kj^jly) 
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Figure  13:  Fdge  stress  intensity  functions  h'a\y)  for  different  degrees  p  at  critical 
edge  labelled  /?i  in  frame  fuselage 
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4  THE  DAMAGE  TOLERANCE  ANALYSIS  OF  COM¬ 
PLEX  THREE-DIMENSIONAL  COMPONENTS 


Fuselages  of  the  type  shown  in  figure  10  often  are  designed  according  to  the  design 
code  [1].  This  code  stipulates  that  pre-existing  flaws  in  form  of  close-fitting  cracks 
of  prescribed  shape  and  size  and  with  most  unfavorable  location  and  orientation 
should  be  considered  in  the  design  process.  We  consider  as  an  example  the  damage 
tolerance  analysis  in  case  of  a  flaw  located  at  edge  B. 

By  the  principles  of  section  3  the  maximum  stresses  are  determined.  The  most 
unfavorable  initial  crack  orientation  for  edge  B  is  such  that  the  plane  of  the  crack 
is  perpendicular  to  the  edge.  Based  on  our  knowledge  of  stresses  wTen  no  flaws  are 
present  (figure  12)  we  select  as  the  most  critical  flaw  location  the  point  of  maximum 
stress  at  edge  B. 

The  calculated  maximum  stress  CT]  (MPa)  at  edge  B  can  be  approximated  with. 


cTi  =  450  —  50d  (16) 

where  d  <  '3rnm  is  the  distance  (mm)  perpendicular  to  the  hole  boundary  (hole 
radius  is  20mni).  The  stress  gradient  over  the  web  thickness  is  small  at  the  point  of 
maximum  stress. 

The  crack  propagation  rate  (m/load  cycle)  of  a  point  at  the  crack  front  is  gener¬ 
ally  assumed  to  be  a  function  of  the  local  value  of  the  mode  I  stress  intensity  func¬ 
tion  The  growth  of  the  flaw  is  obtained  from  integration  of  crack-propagation 
"laws"  to  give  the  life  of  the  component  [2j.  In  critical  cases  the  designer  obviously 
has  to  have  an  effective  tool  for  accurate,  efficient  and  reliable  determination  of  the 
mode  1  edge  intensity  function  for  different  crack  dimensions. 


4.1  Mesh  design 

Figure  15  shows  a  detail  of  the  input  data  mesh  used  to  calculate  the  stress  intensity 
function.s.  The  crack  front  is  divided  in  three  sub-intervals.  The  mesh  generator  is 
used  to  generate  four  layers  of  elements  close  to  the  crack  front.  Several  different 
planar  cracks  with  elliptical  crack  fronts  (semi-axes  a  and  h)  have  been  analyzed. 
The  smallest  crack  size  is  a  =  6  =  1.27mm  [1]  and  the  largest  is  a  =  3.81mm,  b  = 
6.35mm.  Different  FE- models  corresponding  to  different  crack-sizes  are  obtained  by 
simply  ’’stretching”  the  local  mesh  having  a  fixed  number  of  elements. 


4.2  Calculated  edge  intensity  functions 
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Domain 

x' 

y' 

p  =  6 

p  =  5 

p  =  4 

p  =  3 

P  =  2 

1 

1.270 

0.000 

592 

596 

594 

570 

527 

1 

1.26S 

0.064 

597 

594 

589 

567 

524 

1 

1.2G4 

0.128 

592 

589 

582 

559 

520 

1 

1.255 

0.192 

584 

582 

574 

549 

514 

1 

1.244 

0.255 

579 

587 

567 

542 

507 

1 

1.230 

0.318 

570 

577 

559 

541 

498 

2 

1.230 

0.318 

574 

572 

562 

530 

509 

2 

1.135 

0.570 

564 

562 

558 

540 

500 

2 

0.9S8 

0.797 

565 

564 

560 

538 

494 

2 

0.797 

0.988 

572 

569 

561 

533 

490 

2 

0.570 

1.135 

585 

580 

566 

537 

489 

2 

0.318 

1.230 

612 

607 

594 

559 

490 

3 

0.318 

1.230 

608 

612 

589 

556 

494 

3 

0.255 

1.244 

620 

614 

601 

563 

506 

3 

0.192 

1.255 

630 

626 

612 

574 

518 

3 

0.128 

1.264 

641 

636 

623 

586 

529 

3 

0,064 

1 ,268 

650 

646 

634 

598 

541 

0.000 

1.270 

651 

654 

645 

610 

552 

Tahlc  12:  Calculated  A'j’*  at  different  locations  (x',?/')  on  semi-circular  crack  front 
a  —  h  —  \  .'2~jrnn 
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Table  12  shows  the  calculated  (10)  edge  stress  intensity  function  A’)'*  for  a  semi¬ 
circular  crack  of  size  a  =  1 .27r7?m  and  6  =  1.27mm.  VV’e  see  that  for  increa.sing  p.  tljt' 
maxiinuni  value  of  A'j’^  converges  quickly,  smoothly  and  from  below.  Our  prar  a! 
e.xpericiice  is  that  this  is  the  case  for  reasonable  meshes. 

W  c  see  tlist  the  maximum  A’)'^  can  be  computed  virtually  by  guaranteed  accu¬ 
racy  on  models  of  this  (or  higher)  geometrical  complexity. 

For  very  small  crack  sizes  the  stress  intensity  function  A'|''  can  be  roughly  esti¬ 
mated  from  ( 16)  and  weight  function  data  given  in  handbooks  [25].  Such  an  estimate 
is  based  on  the  assumptions  that  crack  dimensions  are  negligible  as  compared  to  the 
hole  radius  and  the  web  thickness.  For  a  =  6  =  1.27mr,i  the  handbook  solution  is 
found  to  o^■erestinlate  the  calculated  maximum  value  (p  =  6)  of  A'j'^  with  about 
ION  . 

.At  the  two  points  where  the  crack  front  intersects  with  the  traction  free  surface. 
A’]'  '  should  be  zero.  This  follows  from  (5)  and  (7)  while  0.5478  and  = 

1,  2.  I’he  calculated  values  of  A’j'*  (table  12)  though  are  not  zero  on  the  surfacel 
Reason  is  that  the  boundary  layers  are  very  local.  For  y  small  y  being  the  distance 
tC)  the  \'crtex,  wt-  have 


i^r(yi  = 


(i)  1 


Polynomials  of  low  order  (e.g.  p  =  6  and  the  three  sub-intervals  used)  simply 
cannot  approximate  the  ste(>[)  gradient  (corresponding  to  the  function  of  the 

edge  intensity  function  close  to  the  two  vertices.  Ii.  ca.se  of  the  folded  plate  analysed 
in  .secti(;n  2  the  (dominating)  gradient  was  much  larger  (compare  figure  6)  and  licnco 
ee.'ier  to  recoNcr.  If  data  for  the  edge  stress  intensity  functions  close  to  the  \'ertev 
at'  of  intde.'t  th-'  \ertex-e(.lge  intensity  factors  can  be  extracted  directl}-.  .A 
ptucedui'e  IS  de>crilied  in  [22  . 

T  iL'iite  Iti  shows  tilt.-  calculated  values  of  A’}**  as  function  of  p  for  three  elhptiea! 
cracl  -  with  (ver_\'  rouchly)  uniform  stress  intensities  along  the  crack  fronts.  I  he 
maximum  values  of  the  calculated  edge  intensity  functions  A’}]'  and  I\\]]  fut  (j  =■ 
b  =  l.dTiiiiri  are  r(’!ati\'e!y  small  (about  25  and  10  units  respectively). 
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5  COMPUTATIONAL  ASPECTS  OF  THE  h-  p  VER¬ 
SION  OF  FEM 

The  goal  of  the  h—p  version  computations  presented  is  to  develop  a  tool  for  rtliahh 
and  efficient  determination  of  any  engineering  data  of  interest  in  complex  real-life 
components. 

In  the  present  work  this  is  obtained  by  using  computational  procedures  which: 

•  give  exponential  rates  of  convergence  in  the  energy  norm  for  accuracies  of 
interest 

•  use  proper  postprocessing  techniques,  which  converge  as  the  energy,  lo  calcu¬ 
late  design  parameters  of  interest 

In  previous  sections  we  have  addressed  these  questions  on  special  examples.  In 
this  sectioi^  we  will  discuss  computational  aspects  of  the  h  —  p  version  of  FEM. 
Instead  of  analyzing  theoretically  these  questions  we  will  report  the  basic  timing 
data  for  the  computations  of  the  fuselage  parts. 

Tlie  h  -  p  %-ersion  when  the  mesh  is  properly  designed  can  be  understood  as 
the  p- version  and  implemented  with  uniform  increase  of  the  degree  p  or  adaptive 
increase  of  p.  The  proper  mesh  design  for  the  present  approach  is  obtained  using 
a  coarse  user-prepared  mesh  combined  with  mesh  generation  at  edges  and  vertices 
where  solution  exhibit  boundary  layers. 

The  main  parts  of  the  analysis  are 

a)  construction  of  the  mesh 

b)  computation  of  local  and  global  stiffness  matrices 

c)  solving  linear  equations 

d)  postprocessing  techniques 

If  the  polynomial  orders  p  of  the  basis  functions  are  assigned  in  an  adaptive  mode 
there  also  is. 

e)  compulation  of  error  indicators 

Construction  of  the  mesh.  This  part  of  the  analysis  is,  in  practical  cases, 
the  most  expensive  one.  The  steady  progress  of  CAD-techniques,  improved  mesh 
generators  [26]  and  interfacing  programs  tend  to  reduce  the  cost,  though  3D  analysis 
of  a  complex  component  is  still  a  major  task  today. 
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5.1  Computation  of  element  stiffness  matrices 

Three-dimensional  p-version  elements  have  0{p^)  degrees  of  freedom  [27]  and  con¬ 
sequently  the  element  stiffness  matrices  have  0(p®)  entries.  If  a  standard  Gaussian 
integration  scheme  with  0{p^)  points  is  used,  the  computation  of  element  stiffness 
matrices  require  0{p^)  operations.  With  such  an  integration  scheme  the  relative 
cost  of  computation  of  element  stiffness  matrices  will  be  very  significant. 

A  fast  method  requiring  only  0(p®)  operations  is  used  in  STRIPE.  The  method 
applies  for  3D  finite  elements  of  arbitrary  shape,  arbitrary  (practically)  Gaussian 
integration  order  used,  and  for  linear  and  nonlinear  analysis.  The  reduction  in  the 
number  of  operations  is  due  to  a  decomposition  of  the  integration  scheme.  First 
0[p^)  basic  3D  integrals  are  determined  using  standard  Gaussian  integration.  The 
C7(p®)  entries  in  the  stiffness  matrix  are  obtained  as  a  linear  combination  of  the  basic 
integrals. 

For  3D  brick  elements,  the  speedup  in  Cpu-time  obtained  as  compared  tc  the 
standard  scheme  [27]  requiring  0{p^)  operations,  is  a  factor  6  to  100  for  p=4  to  10 
(if  (p  -f-  2)^  Gaussian  integration  points  are  used).  With  this  integration  scheme  the 
CPU-times  for  computation  of  element  stiffness  matrices  were  found  to  be  negligible 
(T-2  9c  of  total  CPU-time)  in  all  cases  with  uniform  p  studied  here. 


5.2  Solving  linear  equations 

Direct  methods  (as  opposed  to  iterative  methods)  for  solution  of  the  resulting  set 
of  linear  equations  are  computationally  most  efficient  (in  terms  of  Cpu-time  )  for 
small  and  medium  large  problems  (<  10^  degrees  of  freedom,  say).  Reasons  are  that 
in  most  cases  several  loading  cases  (right  hand  sides)  are  of  importance.  For  thin 
shell  structures  the  performance  of  iterative  solvers  available  are  unsatisfactory.  The 
direct  solver  used  in  STRIPE  is  an  out-of-core  envelope,  with  domain  decomposition 
ordering. 

For  practical  problems  having  of  the  order  10^  degrees  of  freedom  about  907c  of 
the  total  CPU-time  is  spent  in  the  (direct)  equation  solution  phase.  This  is  a  similar 
ratio  cis  for  the  /i- version  of  FE.M. 

Table  13  summarizes  data  for  the  computational  resources  needed  to  analyze 
the  fuselage.  CPU-times  given  include,  in  cases  when  a  direct  solver  were  used,  the 
analyses  for  p  =  2,  3, 4  . . .  Pmar  ■ 

Data  given  in  table  13  apply  for  one  loading  case.  Increase  in  CPU-time  is  a  few 
percent  per  loading  case  for  the  direct  solver  and  about  50%  per  loading  case  for  the 
iterative  solver.  For  damage  tolerance  assessment,  several  (about  6-10)  crack  sizes 
need  be  analyzed.  For  the  case  with  the  crack  at  edge  B  solutions  for  ten  different 
crack  sizes  can  be  obtained  with  less  than  a  three-fold  increase  in  CPU-time.  The 
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at  edges 
figure  13 
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at  fillets 
figure  14 

Crack  at 
edge  B 
figure  15 

Entire 
fuselage 
figure  10a 

Half 
Fuselage 
figure  10a 

Pmax 

5 

6 

6 

7 

8 

No  of  Elements 

778 

778 

610 

5510 

2755 

.V  (kdofs) 
CPU-times  (hrs) 

66 

103 

86 

1268 

853 

.CRAY-XMP/IS 

2 

5 

2 

- 

- 

.IBM  6000/550 

- 

- 

5 

36 

20 

Disc  (Gbyte) 

2 

4 

2 

4 

2 

Table  13:  Model  dimensions  and  computational  resources  needed  to  analyze  fuselage 
frame.  Direct  equation  solver  was  used  for  models  consisting  of  less  than  1000  finite 
elements 

saving  in  Cpu-time  is  obtained  by  using  a  substructuring  technique. 

CPU-times  reported  in  table  13  are  of  the  order  a  few  hours  for  problems  having 
<  10^  degrees  of  freedom.  The  error  in  pointwise  stresses,  stress  intensity  factors 
etc.  though  are  exceptionally  small  for  the  pmaj.-lpvels  reported.  From  a  practical 
point  of  view  a  lower  pmaz  could  have  been  used.  The  computational  cost  is  small 
compared  to  the  cost  of  preparing  the  models. 


5.3  Super  computer  performance  and  large  scale  problems 


The  performance  of  the  h  —p  version  will  strongly  depend  on  the  computer  archi¬ 
tecture  used.  Data  reported  in  table  13  apply  to  a  single-processor  supercomputer 
and  a  pov.erful  work  station.  The  h  -  p  version  leads  to  dense  systems  of  linear 
equations  and  a  very  efficient  utilization  of  the  processors.  For  the  work  station 
IBM  COOO/RS-550  we  obtain  speeds  of  50-60  MFLOPS  independent  of  problem  size. 
For  vectorprocessors  (as  the  CRA'^’-XMP)  the  speed  increases  with  the  problem 
size.  For  the  problem  having  103000  degrees  of  freedom  the  average  (for  p=2.3.4 
and  5)  speed  is  135  MFLOPS  during  equation  solution  on  a  single  processor  CR.A5’- 
XMP/IS.  On  a  4-processor  CR.AY-XMP  a  speedup  factor  of  3.3  has  been  measured. 
Note  that  these  performance  figures  apply  to  cases  with  unstructured  meshes  leading 
to  exponential  convergence  rales. 

For  large  problems  having  of  the  order  10®  degrees  of  freedom  the  disc  storage 
needed  become.s  prohibitive  if  direct  solvers  are  used.  The  preconditioned  con¬ 
jugate  gradient  (PCG)  method  is  perhaps  the  most  efficient  iterative  method  for 
solving  structural  mechanics  problems.  The  p-version  provides  in  a  natural  way, 
a  preconditioner  (in  principle  the  p  =  1  approximation)  for  the  conjugate  gradi- 
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ent  method.  Combined  with  domain  decomposition  techniques,  powerful  iterative 
solvers  for  large  scale  analysis  may  be  developed.  The  iterative  solver  [28],  [29]  was 
used  in  the  present  analyses. 

Table  13  briefly  summarizes  computational  resources  needed  to  analyse  the  com¬ 
plete  fuselage.  Obviously  very  complex  problems  may  be  reliably  solved  on  powerful 
work  stations. 

State  of  the  art  parallel  computers  of  today  may  have  hundreds  of  processors 
where  each  individual  processor  has  similar  performance  as  the  single  processor 
work  station  used  here.  On  such  computers  reliable  analyses  of  entire  aircrafts 
might  become  feasible.  Models  of  this  size  and  complexity  are  not  needed  for  reliable 
computation  of  local  stresses  etc.  though  the  concept  of  using  one  large  model  might 
offer  many  advantages  in  an  industrial  environment. 

Postprocessing  techniques.  The  additional  cost  for  computation  of  stress 
intensity  factors  etc.  are  well  below  1%  (per  loading  case)  of  the  time  needed  to 
analyze  the  fuselage  parts. 


5.4  Computation  of  error  indicators 

Results  above  are  all  for  the  case  when  the  polynomial  order  p  is  uniformly  increased 
in  all  elements.  The  used  elemental  shape  functions  are  of  hierarchical  type  (i.e.  of 
nodal,  edge,  face  and  internal  shape  functions  [27])  which  allow  the  flexible  change 
of  the  degrees  of  elements  and  keep  the  conformity.  Thus,  the  adaptive  selection  of 
element  shape  functions  is  possible.  Here  we  report  some  results  when  the  adap¬ 
tive  scheme  is  based  on  the  energy  norm  and  the  local  elemental  error  indicators 
characterize  the  contribution  of  various  shape  functions  to  the  energy  error.  The 
notation  of  the  degree  Pmai  now  means  the  maximal  degree  of  any  of  the  used  shape 
functions. 

Tlie  cost  for  calculating  the  error  indicators  is  significant.  If  the  element  stiffness 
matrices  are  available,  the  error  indicators  can  be  calculated  with  a  low  cost.  By 
using  the  fast  method  for  computation  of  stiffness  matrices  the  error  indicators  can 
be  calculated  with  a  low  cost. 

As  an  example  we  consider  the  analysis  of  the  stresses  at  the  fillets  with  the  778 
element  mesh.  Solutions  for  uniform  p=6  and  self-adaptive  solutions  to  pmoj-=9  are 
derived.  The  adaptive  scheme  used  is  described  in  [30], 

The  time  needed  for  computation  of  error  indicators  and  (all)  element  stiffness 
data  for  the  self-adaptive  Pmor=9  solution  is  10%  of  the  total  CPU-time  needed. 
Table  14  gives  the  energies  of  the  solutions. 

W'e  note  that  the  energy  convergence  with  the  self-adaptive  solution  scheme  is 
fast  in  this  case.  The  adaptive  scheme  gives,  for  example,  with  25000  degrees  of 
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p 

N 

Mhm 

N 

2 

12480 

1.7828 

8247 

1.7798 

3 

22068 

1.8112 

9745 

1.8023 

4 

39675 

1.8281 

12886 

1.8204 

5 

65676 

1.8356 

16107 

1.8287 

6 

102894 

1.8394 

20133 

1.8348 

7 

25166 

1.8404 

8 

31457 

1.8461 

9 

39321 

1.8519 

Table  14:  Energies  in  solutions  for  778-element  model  analyzed  using  uniform 
p-extensions  and  p-adaptive  extensions.  Data  given  are  scaled  with  a  factor  10“® 

freedom  an  energy  error  which  is  smaller  than  that  obtained  for  the  uniform  p  =  6 
solution  which  has  103000  degrees  of  freedom.  Further,  the  self-adaptive  scheme 
requires  approximately  less  than  half  the  disk  space  needed  for  a  uniform  scheme. 

In  general  the  error  in  energy  norm  is  in  itself  of  less  interest  but  is  a  good  too! 
for  degree  selection  p  if  combined  with  the  accuracy  checks  for  data  of  interest.  The 
CPl’-time  needed  is  a  more  relevant  measure  of  the  effectiveness  than  the  number 
cf  degrees  of  freedom  used. 

In  figure  17  we  show  the  calculated  von  Mises  stress  at  the  two  most  critical 
points  as  function  of  the  p-level  for  the  cases  with  uniform  and  self-adaptive  solution 
schemes.  The  figure  shows  that  considerable  savings  in  CPU-time  might  be  obtained 
by  using  th"  self-adaptive  solution  scheme. 
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6  CONCLUDING  REMARKS 


W’e  have  shown  effective  approaches  for  reliable  analysis  of  complex  three-dimensional 
aircraft  components  with  virtually  guaranteed  accuracy.  The  main  tool  was  the  h  -p 
version  of  the  finite  element  method  combined  with  new  theoretical  approaches. 

We  see  that 

a)  Reliable  and  accurate  accuracy  assessment  for  any  computational  data  of  in 
terest,  for  example  the  stress  intensity  factors  and  functions,  stresses  etc.  is 
available. 

b)  The  design  of  the  mesh  is  very  different  than  that  for  the  /i-version.  Tin 
mesh  has  to  be  strongly  refined.  Refinement  of  a  geometrical  mesh  is  rec¬ 
ommended.  The  construction  of  the  mesh  relatively  simple  when  some  basic 
rules  are  obeyed.  These  rules  can  relatively  easy  be  implemented  in  a  3D  mesh 
generator. 

c)  The  used  degrees  p  of  the  elements  are  in  the  range  4-8  and  the  rate  of  conver¬ 
gence  in  any  data  of  interest  is  exponential  in  the  range  of  practical  accuracy, 

d)  Computation  of  local  element  stiffness  matrices  could  be  costly  if  not  properl> 
implemented.  If  this  computation  is  properly  implemented  experience  show 
that  this  part  of  the  aztalysis  requires  a  small  percentage  of  the  entire  Cfcii- 
lime. 

e)  .Mthough  the  global  stiffness  matrix  is  much  less  sparse  than  the  matrix  it; 
the  classical  h-version.  experience  shows  that  the  computational  cost  for  an 
achieved  accuracy  is  still  very  favourable  compared  with  the  classical  h-version. 
When  a  rdiahU  and  accurate  accuracy  assessment  of  any  data  of  interest  is 
required  or  when  high  accuracy  is  needed  the  h  —  p  version  is  the  only  practical 
alternative. 

fj  .A  direct  equation  solution  technique  is  preferable  for  problems  having  up  to 
50000  to  100000  degrees  of  freedom  especially  when  several  right  hand  sides 
(loads)  are  present.  Powerful  iterative  solvers  which  are  based  on  the  p- version 
techniques  and  the  PCG-method  may  be  used  to  efficiently  solve  complex 
problems  having  several  million  degrees  of  freedom  on  todays  supercomputers. 

g)  New  methods  for  the  computation  of  the  stress  intensity  functions  and  fac¬ 
tors  and  stress  concentration  factors  were  proposed,  thoroughly  analysed  and 
implemented.  The  rate  of  convergence  is  as  the  rate  of  the  energy,  i.e.  as 
the  square  of  the  convergence  rate  of  the  error  measured  in  energy  norm.  In 
the  following  papers  [5], [6], [7],  [8]. [9], [10]  we  elaborate  in  detail  on  these  ap¬ 
proaches. 


FFA  TN  1992-17 


■IS 

7  REFERENCES 

[1]  US  Air  Force.  Airplane  damage  tolerance  requirements  MIL-A-83444  (USAF). 
1974. 

[2]  S.  Suresh.  Fatigue  of  Materials.  Cambridge  University  Press,  1991. 

[3]  C.  Mestenyi,  A.  Miller,  and  W.  Scymezak.  FEARS;  Details  of  mathemathical 
formulation  UNIVAC  1100.  Technical  report  BN-994,  University  of  Maryland, 
1982. 

[4]  R.  E.  Bank.  PLTMG;  A  software  package  for  solving  elliptical  partial  differen¬ 
tial  equations.  Users  Guide  6.0,  SI.AM  Philadelphia,  1990. 

[5]  B.  .^ndersson,  I.  Babuska,  and  T.  Petersdorff.  Computation  of  vertex  singular¬ 
ities  and  intensity  factors  for  Laplace  equation,  (to  appear),  1992. 

[Gj  I.  Babuska,  T.  Petersdorff,  and  B.  Andersson.  Numerical  treatment  of  ver¬ 
tex  singularities  and  intensity  factors  for  mixed  boundary  value  problems  for 
Laplace  equation,  (to  appear).  1992. 

[7]  I.  Babuska,  T.  Petersdorff,  and  B.  Andersson.  Numerical  treatment  of  vertex 
singularities  and  intensity  factors  for  elasticity  equations,  (to  appear),  1993. 

fSj  B.  .Andersson,  1.  Babuska,  and  T.  Petersdorff.  Computation  of  vertex  singu¬ 
larities  and  intensity  factors  for  equations  of  three-dimensional  elasticity,  (to 
appear).  1992. 

[9]  I.  Babuska,  T.  Petersdorff,  and  B.  Andersson.  Computation  of  edge  stress 
intensity  functions  for  Laplace  equation,  (to  appear),  1993. 

[lOj  I.  Babuska,  T.  Petersdorff.  and  B.  Andersson.  Computation  of  edge  stress 
intensity  functions  for  equations  of  three-dimensional  elasticity,  (to  appear). 
1993. 

[11]  K-J.  Bathe,  N.  S.  Lee,  and  M.  L.  Buchalem.  On  the  use  of  hierarchical  mod¬ 
els  in  engineering  analysis.  In  Proceedings  of  the  Workshop  on  Reliability  in 
Computational  Mechanics,  pages  1-1,  1990. 

[12]  1.  Babuska.  The  problem  of  modeling  the  elastomechanics  in  engineering.  Com¬ 
puter  Methods  in  Applied  Mechanics  and  Engineering,  82:155-182,  1990. 

[13]  T.  V  Petersdorff.  Randwertprobleme  der  Elastizitatstheorie  fiir  Polycdersingu- 
laritdten  und  Approximation  mii  Randelementmethoden.  PhD  thesis,  Technical 
University  of  Darmstadt,  Germany,  1989. 

[14]  M.  L.  Williams.  Stress  singularities  resulting  from  various  boundary  conditions 
in  angular  corners  of  plates  in  extension.  J.  Applied  Mechanics,  19;526-528, 
1952. 


FFA  TN  1992-17 


49 


[15]  V.  Z.  Parton  and  P.  ].  Perlin.  Mathemaiical  Methods  of  the  Theory  of  Elasticity. 
Mir,  1984. 

[16]  M.  R.  Dorr.  The  approximation  theory  for  the  p- version  of  the  finite  element 
method.  SIAM  J.  on  Numerical  Analysis,  21:1181-1207,  1984. 

[17]  M.  R.  Dorr.  The  approximation  of  solutions  of  elliptical  boundary  value  prob¬ 
lems  via  the  p-version  of  the  finite  element  method.  SIAM  J.  on  Numerical 
Analysis,  23:58-77,  1986. 

[18]  B.  Guo  and  I.  Babuska.  The  h-p  version  of  the  finite  element  method,  part  1: 
The  basic  approximation  results.  Computational  Mechanics,  1:21^1,  1986. 

[19]  B.  Guo  and  I.  Babuska.  The  h-p  version  of  the  finite  element  method,  part  2: 
General  results  and  applications.  Computational  Mechanics,  1:203-220.  1986. 

[20]  \V.  Gui  ajid  I.  Babuska.  The  h,p  and  h-p  versions  of  the  finite  element  method 
in  1  dimension,  part  1:  The  error  analysis  of  the  p-version.  Numerischt  Math- 
ematik.  49:577-612,  1986. 

[21]  \V.  Gui  and  I.  Babuska.  The  h,p  and  h-p  versions  of  the  finite  element  method 
in  1  dimension,  part  2:  The  error  analysis  of  the  h- and  h-p  version.  Numerische 
Mathematik.  49:613-657.  1986. 

[22]  B.  Andersson,  1.  Babuska.  and  U.  Falk.  Accurate  and  reliable  determination  of 
edge  and  vertex  stress  intensity  factors  in  three-dimensional  elastomechanics. 
In  ICAS-90-4-^-2-  Proceedings  of  1 7th  Congress  of  the  International  Council  of 
the  Aeronautical  Sciences,  pages  1730-1746,  1990. 

[23]  I.  Babuska,  T.  Petersdorff.  and  B.  Andersson.  Asymptotics  of  solutions  of 
boundary  value  problems  at  rounded  corners,  (to  appear),  1993. 

[24]  \V.  J.  Gordon  and  Ch.  A.  Hall.  Transfinite  element  methods:  Blending  function 
interpolation  over  arbitrary  curved  element  domains.  Int.  J.  Num.  Meth.  Eng.. 
21:109-129.  1973. 

[25]  Murakami.  Stress  Intensity  Factors  Handbook.  Pergamon  Press,  1982. 

[26]  M.  S.  Shephard  and  M.  K.  Georges.  Automatic  three-dimensional  mesh  gen¬ 
eration  by  the  finite  octree  technique.  Int.  J.  Num.  Meth.  Eng.,  32:709-739, 

1991. 

[27]  B.  Szabo  and  I.  Babuska.  Finite  Element  Analysis.  J  Wiley  SONS,  1991. 

[28]  J.  Mandel.  Iterative  solvers  for  p-version  finite  elements  in  three  dimensions. 
Prepared  for  the  proceedings  of  ICOSAHOM  92,  Montpellier,  France,  June 

1992,  1992. 


50 


FFA  TN  199J-17 


[29]  J.  Mandel.  Adaptive  iterative  solvers  in  finite  elements.  In  Manolis  Pa- 
padrakakis,  editor,  Solving  Large  Scale  Problems  in  Mechanics:  Development 
and  Application  of  Computational  Solution  Methods.  John  Wiley  and  Sons,  Ltd, 
Chichest, England,  1992. 

[30]  B.  Andersson  and  U.  Falk.  Self-adaptive  analysis  of  three-dimensional  struc¬ 
tures  using  a  p- version  of  finite  element  method.  Report  FFA  TN  1987-31,  The 
Aeronautical  Research  Institute  of  Sweden,  1987. 


FFATN  1992-  17 


issuing  organisaiion 

The  Aeronautical  Research  Institute  of 
Sweden 

(FFA)  P.O.  Box  11021 
S-161  11  BROMMA,  Sweden 


Document  No 

FFA  TN  1992-  17 

Dale 

Aug  1992 

Security 

Unclassified 

Copy  No 

Heg.  No 

No  ot  pages 

50 

yioiect  NO 

HU-3146 

HU-3328 

Order/contract 

Sponsoring  agency 

The  Defence  Materiel  Administration 
(FMV) 

US  office  of  Naval  Research, 
grant  N00014-90-J-1030 


flfle 


Reliable  Stress  and  Fracture  Mechanics  Anah  sis  of  Complex  Aircraft  Components 
using  a  h-p  Version  of  FEM _ 


Authoffs) 

V 

Bcirje  Andersson,  Ivo  Babuska, 

Tobias  von  Petersdorff  and  Urban  Falk 


j^p^jroved  by 


Checked  ty 


iXi 

Martin  Svenzon 

Head  Structures  Department 


ABsiracl” 


We  develop  effective  approaches  with  which  complex  three-dimensional  compo¬ 
nents  may  be  analysed  with  a  high  and  virtuadly  guaranteed  accuracy.  The  main 
coiTiputational  tool  is  a  /i  —  p  version  of  FEM  practically  realised  with  the  p- version 
program  STRIPE  having  a  mesh  generator  for  automatic  mesh  refinement  at  edges 
and  vertices.  Use  of  advanced  extraction  methods  and  new  theoretical  approaches 
give  exponential  convergence  rates  for  accuracies  in  all  engineering  data  of  interest. 
New  methods  for  reliable  calculation  of  local  stresses  and  stress  intensity  data  at 
edges  and  vertices  to  be  used  for  fatigue  dimensioning  at  fillets,  damage  tolerance 
assessment  of  three-dimensional  flaws  etc.  are  given.  A  complex  read-life  problem  is 
reliably  analysed  in  order  to  demonstrate  the  practical  usefulness  of  the  procedures 
advocated.  The  technical  details  are  given  in  forthcoming  papers. 
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